#include "mex.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

int Umap(int i, int j , int NIue, int NJue, int NIut){
    int index;
    int block1,is,js;
    
    if ( j <= NJue )
        index = i + (j-1) * NIut;
    else{
        block1 = NIut*(NJue);
        is = i-NIue;
        js = j-NJue;
        index = block1 + is + (NIut-NIue) * (js-1);
        if (is<1 || js<1)
            mexErrMsgTxt("Umap : Access to U(i,j) out of bounds. \n");
    }
    return index-1;
}

int Vmap(int i, int j, int NIve, int NJve, int NIvt){
    int index;
    int block1,is,js;
    
    if ( j <= (NJve) )
        index = i + (j-1) * NIvt;
    else{
        block1 = NIvt*(NJve);
        is = i-NIve;
        js = j-NJve;
        index = block1 + is + (NIvt-NIve) * (js-1);
        if (is<1 || js<1)
            mexErrMsgTxt("Vmap : Access to V(i,j) out of bounds. \n");
    }
    return index-1;
}

/* 
 * [ Ay ] = getAyOpt( i, j, U, V, NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h)
 * 
 * GETAY Returns the value of Ay at the point (i,j+1/2) on the pressure grid.
 * 
 *
 *  By GDA & CAB
 *  Reviewed by GDMT.
 */

void mexFunction( int nlhs, mxArray *plhs[], 
		 int nrhs, const mxArray *prhs[] )
{
   int i, j;
   int NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt;
   
   double *U, *V;
   double h, V_cen, V_right, V_left, V_top, V_bot, U_tr, U_tl, U_br, U_bl;
   
   double *Ay;
   
   const mwSize *dimsU;
   const mwSize *dimsV;
   mwSignedIndex dimU,dimV,outputSize[1]; 
      
   outputSize[0] = 1;
   
   /* Check for proper number of arguments */
   if ( nlhs != 1) 
      mexErrMsgTxt("Wrong number of output arguments.");
   if ( nrhs != 13 )
      mexErrMsgTxt("Wrong number of input arguments.");
   
   /* Create matrix for the return argument */
   dimsU = mxGetDimensions(prhs[2]);
   dimsV = mxGetDimensions(prhs[3]);
              
   /* Assign pointers to each input and output */
   i = *(double *)mxGetPr(prhs[0]); 
   j = *(double *)mxGetPr(prhs[1]);
   U = (double *)mxGetPr(prhs[2]); 
   V = (double *)mxGetPr(prhs[3]);
   NIue = *(double *)mxGetPr(prhs[4]);
   NJue = *(double *)mxGetPr(prhs[5]);
   NIut = *(double *)mxGetPr(prhs[6]);
   NJut = *(double *)mxGetPr(prhs[7]);
   NIve = *(double *)mxGetPr(prhs[8]);
   NJve = *(double *)mxGetPr(prhs[9]);
   NIvt = *(double *)mxGetPr(prhs[10]);
   NJvt = *(double *)mxGetPr(prhs[11]);
   h = *(double *)mxGetPr(prhs[12]);

   dimU = (mwSignedIndex)dimsU[0];
   dimV = (mwSignedIndex)dimsV[0];

   plhs[0] = mxCreateNumericArray(1, outputSize, mxDOUBLE_CLASS, mxREAL);
   Ay = (double *) mxGetPr(plhs[0]);   
/*   printf("Parameters space [%d,%d,%d,%d,%d,%f,%f]\n",i,j,
           NJue,NIue,NJut,h,U[1]);
*/
/*   #pragma omp parallel for private(xt,yt,ang,dAng,frame,maxEnergy) default(shared) num_threads(omp_get_max_threads())
*/
    V_cen  = V[ Vmap(i  ,j, NIve, NJve, NIvt ) ];

    if (i==NIvt)   
        V_right = V_cen; /*Zero gradient condition*/
/*        V_right = - V_cen; /* */
    else 
        V_right = V[ Vmap (i+1,j, NIve, NJve, NIvt ) ];   

    if (j<=NJve)
        if (i==1) 
            V_left = -V_cen; 
        else 
            V_left = V[Vmap (i-1,j, NIve, NJve, NIvt ) ];
    else
        if (i==NIve+1) 
            V_left  = -V_cen; 
        else 
            V_left  = V[Vmap (i-1,j, NIve, NJve, NIvt )];

    V_bot = V[Vmap (i,j-1, NIve, NJve, NIvt)];
    V_top = V[Vmap (i,j+1, NIve, NJve, NIvt)];

    U_tr    = U[Umap ( i+1,j  , NIue, NJue, NIut)];
    U_tl    = U[Umap ( i  ,j  , NIue, NJue, NIut)];
    U_br    = U[Umap ( i+1,j-1, NIue, NJue, NIut)];
    U_bl    = U[Umap ( i  ,j-1, NIue, NJue, NIut)];

    *Ay  = + ( U_tr + U_br ) * ( V_cen + V_right )
        + ( V_cen + V_top ) * ( V_cen + V_top )
        - ( U_tl + U_bl ) * ( V_cen + V_left )
        - ( V_cen + V_bot ) * ( V_cen + V_bot );

    *Ay  = *Ay/(4*h);

    return;
}


